# This file re-runs BG results
# and our replication of BG at agency and county level
# and stores all of the results
library(plm) #for panel data
library(lfe) #for felm function
library(rlang)
library(tidyverse)
library(broom)
library(purrr)
library(texreg)
options(scipen=0, digits=7) #to standardize number of digits when exporting data. all data here is run using these options.

source("lib/model-functions-BG.R") #these have all the different models we need to run

controlled_subset <- FALSE
year_subset <- FALSE
city_outliers_eliminated <- TRUE
crime_problem_obs_dropped <- FALSE

# load data
if(controlled_subset){
  load('data/city-aggregate-controlled.RData')
  load('data/county-aggregate-controlled.RData')
  load('data/BG_milit.RData') # BG data for exact replication
} else {
  load('data/city-aggregate.RData')
  load('data/county-aggregate.RData')
  load('data/BG_milit.RData') # BG data for exact replication  
}

# restrict data
if(year_subset){
  agency <- agency[agency$year >=2010 & agency$year<=2012,]
  county <- county[county$year >=2010 & county$year <=2012,]
  BG_milit <- BG_milit[BG_milit$year>=2010 & BG_milit$year<=2012,]
} 

if(city_outliers_eliminated){
  agency <- agency %>% filter(allcrime_rate < 1000000)
}

#drop problem observations for UCR data
if(crime_problem_obs_dropped){
  agency <- agency %>% select(-allcrime_rate) %>% rename(allcrime_rate=allcrimewithna_rate) #this drops allcrimerates that are NA instead of counting them as zero
  county <- county %>% filter(lessthan100coverage==0) #this drops all county obs that have less than 100% coverage from ORIs in that county-year
} 

# create panel data
agency <- pdata.frame(agency,index=c("g","year"))
county <- pdata.frame(county,index=c("g","year"))

#make it panel and make variables factors
BG_milit$State <- as.factor(BG_milit$State)
BG_milit$fips <- as.factor(BG_milit$fips)

#make data frame panel
BG_milit <- pdata.frame(BG_milit, index = c("fips", "year"))

## OLS results ---------
OLS.totalcrime.agency <- run.OLS.model.agency(allcrime_rate, agency)
OLS.totalcrime.county <- run.OLS.model.county(allcrime_countysumrate, county)
#----------------------------------------------------------#
## Replication Models for BG using our Agency level data-----
#-----------------------------------------------------------#
allcrime.agency <- run.IV.model.agency(allcrime_rate, agency)
homicide.agency <- run.IV.model.agency(Murder_rate, agency)
robbery.agency <- run.IV.model.agency(Robbery_rate, agency)
assault.agency <- run.IV.model.agency(Assault_rate, agency)
burglary.agency <- run.IV.model.agency(Burglary_rate, agency)
larceny.agency <- run.IV.model.agency(Larceny_rate, agency)
vehicle.agency <- run.IV.model.agency(MV_theft_rate, agency)

## run reduced form models for BG at Agency level ----
allcrime.red.agency <- reduced.agency(allcrime_rate, agency)
homicide.red.agency <- reduced.agency(Murder_rate, agency)
robbery.red.agency <- reduced.agency(Robbery_rate, agency)
assault.red.agency <- reduced.agency(Assault_rate, agency)
burglary.red.agency <- reduced.agency(Burglary_rate, agency)
larceny.red.agency <- reduced.agency(Larceny_rate, agency)
vehicle.red.agency <- reduced.agency(MV_theft_rate, agency)

#------------------------------------------------------------#
## Replication Models for BG Using our County level data -----
#------------------------------------------------------------#
allcrime.county <- run.IV.model.county(allcrime_countysumrate, county)
homicide.county <- run.IV.model.county(MURDER_countysumrate, county)
robbery.county <- run.IV.model.county(ROBBERY_countysumrate, county)
assault.county <- run.IV.model.county(AGASSLT_countysumrate, county)
burglary.county <- run.IV.model.county(BURGLRY_countysumrate, county)
larceny.county <- run.IV.model.county(LARCENY_countysumrate, county)
vehicle.county <- run.IV.model.county(MVTHEFT_countysumrate, county)

# reduced form
allcrime.red.county <- reduced.county(allcrime_countysumrate, county)
homicide.red.county <- reduced.county(MURDER_countysumrate, county)
robbery.red.county <- reduced.county(ROBBERY_countysumrate, county)
assault.red.county <- reduced.county(AGASSLT_countysumrate, county)
burglary.red.county <- reduced.county(BURGLRY_countysumrate, county)
larceny.red.county <- reduced.county(LARCENY_countysumrate, county)
vehicle.red.county <- reduced.county(MVTHEFT_countysumrate, county)


#------------------------------------------------------------#
## Direct replication of BG Models -------
#------------------------------------------------------------#
ols.BG.totalcrime <- run.ols.model(crime,BG_milit) #old model

allcrime.BG <- run.IV.model.BG(crime,BG_milit)
homicide.BG <- run.IV.model.BG(murder,BG_milit)
robbery.BG <- run.IV.model.BG(robbery, BG_milit)
assault.BG <- run.IV.model.BG(assault, BG_milit)
burglary.BG <- run.IV.model.BG(burglary, BG_milit)
larceny.BG <- run.IV.model.BG(larceny, BG_milit)
vehicle.BG <- run.IV.model.BG(mvtheft, BG_milit)

#reduced form
allcrime.red.BG <- reduced.BG(crime,BG_milit)
homicide.red.BG <- reduced.BG(murder,BG_milit)
robbery.red.BG <- reduced.BG(robbery, BG_milit)
assault.red.BG <- reduced.BG(assault, BG_milit)
burglary.red.BG <- reduced.BG(burglary, BG_milit)
larceny.red.BG <- reduced.BG(larceny, BG_milit)
vehicle.red.BG <- reduced.BG(mvtheft, BG_milit)

#------------------------------------------------------------#
## First stage results -------
#------------------------------------------------------------#
BG_first <- summary(allcrime.BG$stage1)$coefficients
BG_agency_first <- summary(allcrime.agency$stage1)$coefficients

BG_county_first <-summary(allcrime.county$stage1)$coefficients

# Save all results ------
# Agency and reduced form
BG_agency_list <- list(allcrime.agency, homicide.agency, robbery.agency, assault.agency, burglary.agency, larceny.agency, vehicle.agency)
BG_agency_results <- tibble(agency = map(BG_agency_list, tidy))
BG_agency_red_list <- list(allcrime.red.agency, homicide.red.agency, robbery.red.agency, assault.red.agency, burglary.red.agency, larceny.red.agency, vehicle.red.agency)
BG_agency_red_results <- tibble(red_agency = map(BG_agency_red_list, tidy))
# County and reduced form
BG_county_list <- list(allcrime.county, homicide.county, robbery.county, assault.county, burglary.county, larceny.county, vehicle.county)
BG_county_results <- tibble(map(BG_county_list, tidy))
BG_county_red_list <- list(allcrime.red.county, homicide.red.county, robbery.red.county, assault.red.county, burglary.red.county, larceny.red.county, vehicle.red.county)
BG_county_red_results <- tibble(map(BG_county_red_list, tidy))
# BG replication and reduced form
BG_list <- list(allcrime.BG, homicide.BG, robbery.BG, assault.BG, burglary.BG, larceny.BG, vehicle.BG)
BG_results <- tibble(map(BG_list, tidy))
BG_red_list <- list(allcrime.red.BG, homicide.red.BG, robbery.red.BG, assault.red.BG, burglary.red.BG, larceny.red.BG, vehicle.red.BG)
BG_red_results <- tibble(map(BG_red_list, tidy))
#Include OLS
BG_OLS_results <- tibble(map(list(OLS.totalcrime.agency,OLS.totalcrime.county, ols.BG.totalcrime),tidy))
#First stage
first_stage <- list(BG = BG_first, rep_agency = BG_agency_first, rep_county = BG_county_first)
first_stage_results <- tibble(first_stage = map(first_stage, tidy))

# save
if(controlled_subset==F & year_subset==F & city_outliers_eliminated==T & crime_problem_obs_dropped==F){
  save(BG_agency_results,
       BG_agency_red_results,
       BG_county_results,
       BG_county_red_results,
       BG_results,
       BG_red_results,
       BG_OLS_results,
       first_stage_results,
       file = 'data/results-BG.RData')
}

# Save all results for texreg------
# Agency and reduced form
BG_agency_results <- map(list(allcrime.agency, homicide.agency, robbery.agency, assault.agency, burglary.agency, larceny.agency, vehicle.agency), 
                         texreg::extract)
BG_agency_red_results <- map(list(allcrime.red.agency, homicide.red.agency, robbery.red.agency, assault.red.agency, burglary.red.agency, larceny.red.agency, vehicle.red.agency),
                             texreg::extract)
# County and reduced form
BG_county_results <- map(list(allcrime.county, homicide.county, robbery.county, assault.county, burglary.county, larceny.county, vehicle.county), 
                         texreg::extract)
BG_county_red_results <- map(list(allcrime.red.county, homicide.red.county, robbery.red.county, assault.red.county, burglary.red.county, larceny.red.county, vehicle.red.county), 
                             texreg::extract)
# BG replication and reduced form
BG_results <- map(list(allcrime.BG, homicide.BG, robbery.BG, assault.BG, burglary.BG, larceny.BG, vehicle.BG), 
                  texreg::extract)
BG_red_results <- map(list(allcrime.red.BG, homicide.red.BG, robbery.red.BG, assault.red.BG, burglary.red.BG, larceny.red.BG, vehicle.red.BG), 
                      texreg::extract)
#Include OLS
BG_OLS_results <- map(list(OLS.totalcrime.agency,OLS.totalcrime.county, ols.BG.totalcrime), 
                      texreg::extract)
#First stage
BG_first <- allcrime.BG$stage1
BG_agency_first <- allcrime.agency$stage1
BG_county_first <- allcrime.county$stage1
first_stage_results <- list(BG = BG_first, rep_agency = BG_agency_first, rep_county = BG_county_first)

#save
if(controlled_subset==F & year_subset==F & city_outliers_eliminated==T & crime_problem_obs_dropped==F){
  save(BG_agency_results,
       BG_agency_red_results,
       BG_county_results,
       BG_county_red_results,
       BG_results,
       BG_red_results,
       BG_OLS_results,
       first_stage_results,
       file = 'data/results-BG-for-tables.RData')
}

if(controlled_subset==T){
  save(BG_agency_results,
       BG_agency_red_results,
       BG_county_results,
       BG_county_red_results,
       BG_results,
       BG_red_results,
       BG_OLS_results,
       first_stage_results,
       file = 'data/results-BG-for-tables-controlled.RData')
}

if(year_subset==T){
  save(BG_agency_results,
       BG_agency_red_results,
       BG_county_results,
       BG_county_red_results,
       BG_results,
       BG_red_results,
       BG_OLS_results,
       first_stage_results,
       file = 'data/results-BG-for-tables-year-subset.RData')
}

if(city_outliers_eliminated==F){
  save(BG_agency_results,
       BG_agency_red_results,
       BG_county_results,
       BG_county_red_results,
       BG_results,
       BG_red_results,
       BG_OLS_results,
       first_stage_results,
       file = 'data/results-BG-for-tables-with-outliers.RData')
}

if(crime_problem_obs_dropped==T){
  save(BG_agency_results,
       BG_agency_red_results,
       BG_county_results,
       BG_county_red_results,
       BG_results,
       BG_red_results,
       BG_OLS_results,
       first_stage_results,
       file = 'data/results-BG-for-tables-crime-obs-dropped.RData')
}


